A machine learning approach in the non-invasive prediction of intracranial pressure using Modified Photoplethysmography

The ideal Intracranial pressure (ICP) estimation method should be accurate, reliable, cost-effective, compact, and associated with minimal morbidity/mortality. To this end several described non-invasive methods in ICP estimation have yielded promising results, however the reliability of these techniques have yet to supersede invasive methods of ICP measurement. Over several publications, we described a novel imaging method of Modified Photoplethysmography in the evaluation of the retinal vascular pulse parameters decomposed in the Fourier domain, which enables computationally efficient information filtering of the retinal vascular pulse wave. We applied this method in a population of 21 subjects undergoing lumbar puncture manometry. A regression model was derived by applying an Extreme Gradient Boost (XGB) machine learning algorithm using retinal vascular pulse harmonic regression waveform amplitude (HRWa), first and second harmonic cosine and sine coefficients (an1,2, bn1,2) among other features. Gain and SHapley Additive exPlanation (SHAP) values ranked feature importance in the model. Agreement between the predicted ICP mean, median and peak density with measured ICP was assessed using Bland-Altman bias±standard error. Feature gain of intraocular pressure (IOPi) (arterial = 0.6092, venous = 0.5476), and of the Fourier coefficients, an1 (arterial = 0.1000, venous = 0.1024) ranked highest in the XGB model for both vascular systems. The arterial model SHAP values demonstrated the importance of the laterality of the tested eye (1.2477), which was less prominent in the venous model (0.8710). External validation was achieved using seven hold-out test cases, where the median venous predicted ICP showed better agreement with measured ICP. Although the Bland-Altman bias from the venous model (0.034±1.8013 cm water (p<0.99)) was lower compared to that of the arterial model (0.139±1.6545 cm water (p<0.94)), the arterial model provided a potential avenue for internal validation of the prediction. This approach can potentially be integrated into a neurological clinical decision algorithm to evaluate the indication for lumbar puncture.


Introduction
Investigation of the physiological properties of the human cerebrospinal fluid was first described in 1891 when Quinke published his studies on the diagnostic and therapeutic applications of lumbar puncture. He standardized the technique of intracranial pressure (ICP) measurement by connecting the lumbar puncture needle with a fine glass pipette in which cerebrospinal fluid was allowed to rise, a method still currently applied [1]. Whereas data from a recent large international, multi-center study on lumbar puncture feasibility that included 3,868 patients confirmed the procedures' safety [2], complications from lumbar puncture including iatrogenic meningitis, hemorrhage, post-lumbar puncture headache are well recognized [3][4][5][6]. In addition to the risks stated above, continuous ICP monitoring includes ventricular catheter-related problems such as cerebral infections, catheter occlusion, or malposition. [7][8][9]. Financial costs incurred at a single tertiary care institution due to hospitalization for post-lumbar puncture complications were estimated at $20,000 USD/year [10]. To mitigate these risks, various modalities of non-invasive ICP estimation have been described. Of these studies ophthalmodynamometry provided the earliest attempts at non-invasive ICP estimation, they were based on earlier observations by Deyl in 1898 who postulated that papilledema was due to central retinal vein compression where it emerged from the optic nerve into the subarachnoid space in the optic nerve's dural sheath [11]. To further prove this mechanism, Cushing and Borley experimentally induced papilledema and described the loss of spontaneous venous pulsation due to collapse of the central retinal vein when ICP was raised [12]. Their findings were supported by a canine model, in which the temporal succession of events demonstrated that loss of spontaneous venous pulsation and venous dilation preceded optic disc enlargement was demonstrated [13]. In 1927 Baurmann provided the first evidence of a strong linear correlation between ICP and retinal venous pulse pressure measured by ophthalmodynamometry [14]. His results were replicated by other investigators in a series of animal studies [15][16][17][18][19]. These findings were also reported in human studies by Firsching [20,21], Motschmann [22] and co-workers. To improve ICP predictive accuracy, Querfurth et al. combined ophthalmodynamometry with simultaneous color doppler of the central retinal and ophthalmic arterial flow velocities. Although the linear correlation between retinal venous pulse pressure and ICP was strong. The combination of ophthalmodynamometry and color doppler parameters was shown to improve the correlation compared to either parameter alone [23]. Intuitively, a strong linear correlation would imply strong predictive accuracy, yet paradoxically these linear models did not achieve a strong predictive power expected from the linear correlation and none of these methods have superseded invasive methods of ICP measurement to date [24][25][26].
The technique and instrumentation for ophthalmodynamometry was developed by Bajardi around 1906 for the indirect estimation of ocular perfusion [27]. Although this device and its several iterations was superseded by carotid doppler, we have over a series of publications [28-34] described a novel imaging method of Modified Photoplethysmography where a combination of ophthalmodynamometry as a means of generating a range of induced intraocular pressures and slit-lamp imaging of the optic nerve enables modelling of the retinal venous and arterial pulse amplitude and timing characteristics. The physiological basis for the choice of a harmonic regression model in image analysis is to emulate the photoplethysmography wave periodic and non-periodic components [35]. The model consists of a linear spline, which represents the mean of the signal and adjusts for inter-frame image displacement. The first two harmonics of the Fourier series are fitted to the periodic component and a first-order autoregressive error component accounts for the error process. Frequency domain decomposition is performed using custom software, where heat maps of the retinal vascular pulse amplitude distribution are generated and Fourier coefficients of the first two harmonics are extracted. Additional to modelling the retinal arterial and venous systems separately using this method, uniquely, frequency domain analysis allows computationally efficient information filtering and comparative processing of the retinal vascular pulse characteristics [30,36]. The trade-off in signal resolution between the temporal (Δt) and frequency (Δf) domains is defined by the Heisenberg-Gabor uncertainty principle. It is mathematically expressed as an inequality equation (Δf � Δt � C), where (C) is a constant with a value dependent on how the frequency is measured and conveys the general concept that a non-zero function and its Fourier decomposition cannot be localised to arbitrary precision [37]. Although information between these domains is conserved, unique hemodynamic phenomena can exist solely as a property of the frequency domain such as dispersion (frequency dependent velocity of a wave), impedance (frequency dependent resistance), and the harmonic amplitude distribution. Modified photoplethysmography requires imaging at a slit-lamp, therefore is best applied in an ophthalmology outpatient setting, this excludes the applicability of the method for critically ill patients in the emergency or intensive care setting. A sophisticated handheld system is currently under development by our group, which may make future iterations of the clinical system more versatile for ICP measurement. Although it is possible to perform the test by a single operator, we have used an observer to vocalize the force readings from the ophthalmodynamometer, this allows editing the videos as further detailed in the materials and methods section. Moreover, the technique requires patient cooperation to maintain fixation and remain seated as is required for slit-lamp ophthalmoscopy, this excludes patients with cognitive impairment. A reduced model of the pulsation event is possible from the output of Modified Photoplethysmography it cannot provide information on the hemodynamic pressure-flow wave and although it provides heat maps of the pulse amplitude and timing characteristics, further work is required to analyze vascular geometry specific changes of the pulse characteristics in the Fourier domain. The interaction between ICP, IOP, and the retinal vascular pulse may share characteristics of chaotic systems. In general, systems with at least two of the following properties are considered to be chaotic: bifurcation and period doubling, period three, transitivity and dense orbit, sensitive dependence to initial conditions, and expansivity [38]. Even though Modified Photoplethysmography, unlike other non-invasive methods of ophthalmic ICP estimation, provides continuous vascular observations to potentially resolve this question. The limitation in our work is date is due to the lack of concurrent continuous ICP and IOP measurements during retinal imaging, which particularly the latter, may not be possible out of an experimental / physical model setting. In a recently published study, we described the interaction of the harmonic regression amplitude (HRW a ) distribution in the retinal vascular system with intraocular and intracranial pressure using a linear mixed effects model. This approach enabled the computation of the variance estimated by these variables in linear space. It was demonstrated that linear interactions of IOP, ICP and the retinal vascular pulse accounted for less than 10% of the variance. This poor explanatory power of a linear model precludes it as a predictor of ICP. This was due to a non-constant variance of the error term of the predictors (heteroscedasticity), which indicated that the correlation was linear within individuals but non-linear between individuals [28]. Therefore, addressing this non-linear component is crucial to achieving high ICP predictive accuracy. In this study, we applied an Extreme Gradient Boost (XGB) supervised machine learning decision tree algorithm in noninvasive ICP estimation from features extracted from the retinal arteries and veins separately. This approach addresses the non-linear component in the data structure and model fit was adjusted through tuning of the hyperparameters. Additionally, parallelization and distributed computing allowed for rapid data computation run times. We hypothesize that the combination of Fourier-domain decomposition of the retinal vascular pulse wave parameters and a regression model derived from a decision tree machine learning algorithm could provide an avenue to non-invasive ICP prediction.

Subject recruitment
Twenty-eight participants were recruited prospectively from the Lions Eye Institute over five years (2015-2020) from referrals made to the clinic for ophthalmic assessment before lumbar puncture for suspicion of idiopathic intracranial hypertension. Study approval was obtained from the University of Western Australia Human Ethics Committee adhering to the tenets of the Declaration of Helsinki. Participants were required to have clear ocular media, no prior history of co-existing retina or optic nerve disease, and were needed to be able to cooperate with the imaging protocol. Written consent was obtained from each of the participants. Lumbar puncture was performed in the lateral decubitus position and measured in centimeter (cm) water. The ophthalmic examination consisted of measurement of visual acuity, Goldmann tonometry, slit-lamp examination, color fundus photography, and modified photoplethysmography. The latter test consists of contact lens ophthalmodynamometry, the purpose was to vary induced intraocular pressure (IOP i ), with concomitant video imaging of the optic disc.
An ICP of 25 cm water was considered the upper normal limit [39], this threshold classified the twenty-one patients in the training and test study groups into ten cases in the high intracranial pressure group (ICP h >25cm water) and eight in the normal intracranial pressure group (ICP n �25cm water). Three cases overlapped both groups as a result of interchanging between the ICP n to the ICP h groups over the observation period (Fig 1). A total of 129,600 data points were sampled from the images, 56,932 arterial and 72,668 venous data points ( Fig  2). Three eyes were excluded from the analysis due to poor image quality. Additionally, seven subjects were recruited for model validation. data from these cases were not subjected to the training or testing phases of the analysis.

Image acquisition
Data capture was performed by an operator and an observer in an outpatient setting. While the operator concentrated on the imaging task, the observer vocalized the ophthalmodynamometric force, this was required to create the splice points for video editing. The test requires patient cooperation to sit upright, remain stationary and fixate on a target as required for slitlamp ophthalmoscopy. Details for the image acquisition workflow is published in earlier work [30]. The optic nerve was imaged under a dynamic range of intraocular pressures using a Meditron ophthalmodynamometer (Meditron GmbH, Poststrasse, Völklingen, Germany). This device consists of a sensor ring, which measures the compression force on the eye. The sensor surrounds a central Goldmann three-mirror fundus contact lens. The ophthalmodynamometric force (ODF) displayed as Meditron units (mu), which were then converted to induced intraocular pressure (IOP i ) as described by Morgan et al. [40] using the following formula: Where IOP b is the baseline intraocular pressure (IOP) in millimeters mercury (mmHg). Video of the optic nerve was captured with an imaging slit-lamp (Carl Zeiss, Germany) with a mounted digital camera (Canon 5D Mark III, Japan). Several sequences of at least three cardiac cycles in length were taken, each at a rate of 25 frames/second. When possible recordings were

PLOS ONE
taken from both eyes. A range of induced intraocular pressure values was between 7-73 mmHg were obtained from each subject. Videos showing motion artifact, reflection from optical media, or decentration of the optic nerve in the image sequence for less than three consecutive cardiac cycles were rejected from the analysis. A pulse oximeter (Nellcor N65, Covidien, Mansfield, MA) was applied to the right index finger; the audio signal from the pulse oximeter was recorded with the video sequence of the optic nerve. This allowed synchronization of the retinal vascular pulse with the cardiac cycle. Timing of the cardiac cycle was generated from the audio signal from the subject's pulse oximetry recorded on the audio trace of the video segment, which in turn enabled the mathematical analysis of the periodic component from green channel transmittance. A Single high-quality three-cardiac cycle length video recording was extracted from each recording session.

Image analysis
Image processing was done in Adobe Photoshop CS6, Individual image frames were extracted from each video sequence and saved as Tagged Image File Format (TIFF) files. Each of these images was cropped to an array of pixels. All images from three cardiac cycles were analyzed in R statistical package using custom software [41]. Each data point was represented by the mean of the green channel intensity at time measured as a fraction of the cardiac cycle, rather than in seconds. The periodic trend component was modelled separately for the arteries and veins as a harmonic regression waveform expansion:

Machine learning algorithm-Extreme Gradient Boost (XGB)
Extreme gradient boost algorithm is an ensemble machine learning regression method based on decision trees that use the gradient descent architecture to boost weak learners [42]. Boosting builds decision trees sequentially such that each subsequent tree aims to reduce the errors of the previous tree and the residual errors are then updated. R statistical package [41] was used to generate the model for each vascular system independently. Bayesian optimization was used to tune seven of the model hyperparameters aimed at regulating the model fit. Five-fold cross-validation with ten early stopping rounds was applied in this step. Tuned hyperparameters included lambda (λ) L2 regularization term on weights (analogous to Ridge regression), this parameter reduces the influence of outliers by factoring in the denominator of the similarity score (defined as the ratio of the sum of the residuals squared and the number of the residuals). Gamma (γ) is the minimum loss reduction required to make a further partition on a leaf node of the tree. The larger gamma is, the more conservative the algorithm. Eta (η), the learning rate, is a scalar that determines step size in gradient descent and shrinks the weights on each step. Alpha (α) L1 regularization term on weights (analogous to Lasso regression). The maximum depth of the tree (max_depth) is the maximum number of nodes allowed from the root to the farthest leaf of a decision tree, deeper trees can model more complex relationships by adding more nodes. Minimum child weight (min_child_weight) is the minimum number of samples required to create a new node in the tree. Subsample corresponds to the fraction of observations (the rows) to subsample at each step, and nrounds the number of model iterations [42,43]. Hyperparameters values for each model are listed in Table 1. If parameters are not set, default values are chosen by the XGB algorithm.
In our study the training data consisted of nine numerical features (IOP i , HRW a , the cosine and sine coefficients of the first and second harmonic waves (a n1,2 ), b n1,2 ), hemiretinal location of the vessel (superior or inferior retina), and laterality (right or left eye)). The training labels which were the ICP measured by lumbar puncture in the lateral decubitus position in centimeter (cm) water, were ultimately used to generate the ICP predictions. The training set consisted of 80% of randomly selected vascular pulsation points, the test set consisted of the remaining 20% of the data points. Model parameters were assessed using feature importance, which is a ranking score representing the contribution from the selected feature to the model prediction. It is calculated for a single decision tree by the amount that each attribute split point improves the performance measure, weighted by the number of observations for which the node is responsible. There are three methods for measuring feature importance in XGB, frequency (weight), which is the number of times a feature is used to split the data across all trees. Cover is the number of times a feature is used to split the data across all trees weighted by the number of training data points that go through those splits, and gain is the average training loss gained when using a particular feature at a branching point. Gain, therefore, represents the refinement in accuracy brought by a feature to the branches of the decision tree. XGB divides feature importance by default into two clusters, cluster one contains features of the highest importance to the model, and other features are aggregated in cluster two. To identify the main features driving model prediction, SHAP (SHapley Additive exPlanations) values were calculated, this is an additive feature attribution method that provides a quantitative evaluation of the tree ensemble's overall impact in the form of particular feature contributions [44,45]. External validation was achieved using seven hold-out test cases, which were evaluated in a blinded manner. From the hold-out test set mean, median and peak density (defined as the maximum point generated from the distribution of predicted ICPs from each model) of the predicted ICP were compared to the measured ICP. Both Bland-Altman plots and the t-test were used to measure the mean difference between predicted and measured ICP for both the arteries and the veins separately. A flow chart of image processing, analysis, and XGB model application is summarised in Fig 3.

Statistical analysis
The distribution of the HRW a and the majority of the Fourier coefficients was non-normal, therefore the median was used as a measure of central tendency and the interquartile range  3) The Fourier coefficients, harmonic regression amplitude, together with the distance along the vessel, induced intraocular pressure (IQR) was used to assess the dispersion of this measure. Where appropriate the mean and standard deviation were reported. The range, minimum, and maximum of these parameters were also computed. When the Levene test was applied to assess the multifactorial homogeneity of variance of the predictors, heteroscedasticity was demonstrated as the assumption of homogeneity of variance was violated (p<0.0001). The Kruskal-Wallis test was used in the hypothesis test of the differences in the medians, and the paired Wilcoxon test with Bonferroni-Holm correction was used for posthoc analysis. Model fit was assessed using R 2 square, which is a comparison of the residual sum of squares with the total sum of squares. Model Prediction accuracy was estimated by calculating the Mean Squared Error (MSE), defined as the square of the difference between the predicted and actual values of the test set, it assigns more weight to larger errors. Root mean square error (RMSE), which is the standard deviation of the residuals (prediction errors), the higher the number the greater the standard deviation σ of the distribution of errors. MSE and RMSE are used to evaluate the influence of outliers on predictions.
The mean absolute error (MAE) calculated by the magnitude average difference between the predicted and actual values of the test set was also reported.

Descriptive statistics
There were a total of twenty females (95.2%) and one male (4.8%) in the study population. The age demonstrated a bimodal distribution with a mean of 32 years (sd 8.32, range 17-47 years).
In the ICP n group median ICP was 18.50 cm water (range 9.50 to 24, IQR = 6), the corresponding values in ICP h group were 31 cm water (range 25.50 to 68, IQR = 10). Table 2 demonstrates the Fourier wave amplitude descriptive parameters in both study groups. Hypothesis tests within ICP group differences were contrasted by vessel type, statistical significance (p<0.0001) was achieved for all Fourier parameters except the b n2 coefficient in the ICP h group, which demonstrated no statistically significant difference between the retinal arteries and the retinal veins. Both venous and arterial Fourier wave amplitude descriptive parameters stratified by ICP (between-group differences) demonstrated statistical significance (p<0.001) for all parameters except the retinal venous b n1 , and arterial a n2 coefficients. The ICP h group showed a lower median retinal venous (4.743 vs 5.314) and a higher arterial HRW a (4.559 vs 4.139) compared to the ICP n group (Fig 4).

Machine learning model
Model fit and feature importance. The model fit was was comparable for the arterial and venous data as indicated by an R 2 =0.89 and R 2 =0.91 respectively. Other accuracy parameters were similar for the arterial (MSE = 10.99, MAE = 2.03, RSME = 3.32) and the venous model (MSE = 11.85, MAE = 2.11, RSME = 3.44). The venous decision tree was more complex than the arterial, whereas the venous model was composed of a total of 451 nodes, 450 edges, and 35,589 leaves, the arterial model consisted of a total of 137 nodes, 136 edges, and 7,951 leaves. Model complexity and fit of the arterial and venous models are demonstrated in (Figs 5 and 6). The weighted cover in the figures represents the distribution of the average weighted number of residuals clustered in leaves at a certain depth of the decision tree. Global feature importance ranks the nine features of each vascular model by the feature gain, cover and frequency used in the prediction of ICP (Figs 7 and 8). For the arterial model these were IOP i (0.6092), a n1 (0.1000) and HRW a (0.0804), similarly for the venous model IOP i (0.5476) and a n1 (0.1024), dominated the feature importance, however unlike the arterial model HRW a (0.124) showed higher importance compared to a n1 coefficient. When feature frequencies of the arterial and venous models were compared IOP i (0.2547 vs 0.2028), a n1 (0.1357 vs 0.1362) and HRW a (0.1228 vs 0.1291), accounted for approximately 51-47% of each model's feature importance respectively (Table 3).
Model SHAP values. The SHAP summary plot (Figs 9 and 10) combines feature importance with feature effects, therefore it allows to explore interactions between features for the predicted variable. It is important to consider that SHAP values do not identify causality [45,46]. Four properties can be derived from the SHAP summary plot: 1. Feature importance: features are ranked on the y-axis in descending order according to their importance in the prediction of ICP from each vascular model.
2. Impact: SHAP measures the impact of variables taking into account the interaction with other variables of the model. The horizontal location shows whether the effect of that value is associated with a higher or lower prediction. This is accomplished by calculating the importance of a feature by comparing what a model predicts with and without the feature  From both plots, it can be observed that the impact of IOP i on the model prediction was dependent on the value and showed a positive correlation, low IOP i values had a low predictive impact, and higher IOP i values had a positive impact on prediction, particularly for the venous model. This feature also showed higher dispersion than any of the tested features. The hemiretinal location of the tested vessel showed the lowest impact on model predictability. Laterality and distance of the data point along the vessel measured from the center of the optic disc

PLOS ONE
attained more significance compared to the feature importance without variable interactions (Table 3), particularly for the arterial model, where there was a correlation with laterality, right eyes had negative and left eyes a positive impact on prediction, the venous model did not demonstrate this effect. Pulsation values obtained from vascular points in proximity to the optic disc had a less predictive impact and pulsation values from more peripheral locations in the vessel had a higher impact on model prediction (positive correlation).
The arterial model showed no correlation with HRW a values, furthermore, this feature was less significant when interactions were considered. In contrast, the venous model

PLOS ONE
demonstrated a negative correlation, and it retained its significance both with and without variable interactions ( Table 3).
The a n1 Fourier coefficient had the highest feature importance of all coefficients in both models, and other coefficients had a low rank. The correlation of the coefficients was different for each model. Whereas the arterial model demonstrated that the a n1,2 , and b n1,2 showed low positive correlation on model predictability, the venous model on the other hand showed that a n1,2 and b n1,2 had opposite correlation on predictability, where a n1 , b n1,2 values had a negative correlation and a n2 had a positive correlation.

PLOS ONE
Mean SHAP values are listed in Table 4, where IOP i demonstrates the highest mean SHAP values in both vascular models (arterial = 5.3884 and venous = 5.6375), this was approximately four times the mean value of a n1 (arterial = 1.4689 and venous = 1.3856) and the others among the three most significant features (arterial laterality = 1.2477, venous HRW a =1.7024).

Model external validation
A group of seven cases was used as the holdout test set, in this group all cases were females. Four cases (57.1%) had an ICP > 25cm water. The median ICP in this group was 29.5cm water (range 26 to 32, IQR = 5.25). The remaining three cases (42.9%) had a median ICP of 20cm

PLOS ONE
water (range 17 to 22, IQR = 2.5). Table 5 summarises the mean, median, and peak density of the predicted ICP from each vascular model. When the mean, median, and peak density for predicted ICP from the arterial and venous models are compared using the t-test and Bland-Altman bias statistic (Table 6), it is clear that the predicted ICP estimated from the venous median had the best agreement with measured ICP as indicated by the lowest Bland-Altman bias. A comparison of measured and median estimated ICP is demonstrated graphically in Bland-Altman plots (Figs 11 and 12).

Discussion
Extreme Gradient Boost demonstrated favorable accuracy in the non-invasive prediction of ICP applied to Modified Photoplethysmography data. Quantitative interference due to optical reflections, shadowing, blink, and motion artifact secondary to saccadic movements render ophthalmic imaging artifact prone. Modified Photoplethysmography controls these multiple sources of interference. The Goldmann contact lens used for optic nerve observation and imaging eliminates blinking and reduces motion artifact providing optical continuity, field stability, and reducing information degradation. Induced intraocular pressure generates a range vascular pulse amplitude responses, therefore, enabling comparative analysis of the pulse wave under a range of transmural pressures. Moreover, retinal vascular pulse wave decomposition in the Fourier domain allows for computationally efficient information filtering. The harmonic regression approach applied in image analysis not only adjusts for motion artifacts through its linear spline; it applies a statistical approach to evaluate the fit of the Fourier harmonics to the non-periodic component of the vascular pulse. Hence it facilitates the decision of rejecting an  Over the last two decades, OCT has been central to the diagnosis and management of optic nerve disorders, the earliest report was by Borchert et al. who patented a method to estimate ICP using OCT measurements of RNFL thickness; however, the authors do not provide the correlations between these variables necessary generate a prediction [55]. A multitude of parameters have been evaluated for potential estimation of ICP [56-60]. In a multicenter study, Vijay et al. reported that optic nerve head central thickness was found to be the most closely associated parameter with ICP (R = 0.60-0.73) among a variety of macular and optic nerve protocols [61]. However, the association between the many OCT parameters, papilledema severity, and ICP is complex and remains undefined. Moreover, due to the need for subject cooperation, this test cannot be applied to patients with severe neurological disorders. To address this limitation, Andersen et al. used the retinal arterio-venous ratio as a biomarker of elevated ICP. Images were recorded using an Epicam portable camera, this method achieved a 94% (85-98%) sensitivity and 50% (34-66%) specificity in detecting patients with an The measured ICP was compared against the mean and peak density of the estimate from both vascular models. The venous median demonstrates the highest agreement with measured intracranial pressure. sd = standard deviation, se = standard error. Bland-Altman plots are displayed in (Figs 11 and 12). https://doi.org/10.1371/journal.pone.0275417.t006 ICP � 20mmHg. Indicating that although there was a 94% probability of correctly identifying individuals with ICP � 20 mmHg, this was mitigated by the 50% probability of misclassifying healthy individuals [62], thereby limiting the practical applicability of this approach. Radiological imaging of retinal vascular parameters has been explored as a substitute for clinical methods especially since practicality demands the ability to predict ICP with minimal patient cooperation in the ICU setting. In a prospective case-control study, Jeub et al.

PLOS ONE
threshold value of 11.0 cm/s, the peak systolic velocity predicted pathological ICP levels with a 70% sensitivity and 69% specificity [63]. Using spectral Doppler imaging, Miller et al.

PLOS ONE
this method was reported to have better diagnostic reliability for detecting elevated ICP compared to optic nerve sheath diameter measurement [66][67][68][69]. An independent clinical validation study determined that this technique had a fair agreement to ICP measured using lumbar puncture [70]. The theoretical basis establishing a correlation between optic nerve sheath diameter and ICP was first suggested in 1968 by Hayreh et al. in a primate model [71]. Subsequent studies have imaged this parameter using U/S, CT, and MRI making it a favorable option for subjects with significant neurological impairment. In the largest study to date, Rajajee et al. found the optimal optic nerve sheath diameter for detection of ICP > 20 mmHg was >0.48 cm as measured by U/S, where a sensitivity of 96% and specificity of 94% were achieved [72]. Most studies indicate an optic nerve sheath diameter of >5 mm as the threshold for determining elevated ICP [73][74][75]. Measurement of the optic nerve sheath diameter using CT or MRI can overcome this limitation. High agreement and reproducibility have been reported between these imaging modalities [76][77][78]. However, radiological methods have high operational costs, and image interpretability is operator dependent particularly in the case of U/S. Some disorders such as subarachnoid hemorrhage may not be suited for this modality [79]. Moreover, the range of pressure correlation with optic nerve sheath diameter can be narrow, the iCOP study reported a favorable correlation between optic nerve sheath diameter and ICP between 3.7 mm Hg and 26.5 mm Hg [80]. The precision and accuracy of MRI measurements of optic nerve sheath diameter are yet to be defined, as well as the optimal measurement technique and the influence of the time course of ICP fluctuations on changes in optic nerve sheath dimension [81].
Flash Visual Evoked Potentials (FVEPs) can demonstrate the integrity of the visual pathway from the retina to the occipital cortex. A longer FVEP wave crest latency, a reduction in amplitude, and an increase in wave width have been observed with elevated ICP. These findings were reported from two early studies where a relatively strong linear relationship (R 2 � 0.7) between FVEP N2 wave latency and ICP was observed. Using this method high correlation was particularly demonstrated at ICP levels >300 mm water [82,83]. These findings were replicated using either a single or combined modality with transcranial doppler [84,85]. There are significant limitations with FEVPs, a significant level of expertise is required to administer the test, it is unsuitable in patients with bifrontal lobe pathology, retinal damage, or optic neuropathy [84]. Moreover, factors such as blood glucose concentration, the patient's nerve conduction rate, and electrolytes levels can result in high variance in the FVEP waveform properties [86].
Early studies which applied ophthalmodynamometry to this research question correlated retinal venous pulse pressure as a single parameter with ICP [20][21][22]. The wide range of linear correlations (R = 0.69-0.968) did not, however, yield a clear conclusion regarding the applicability of their findings. Further attempts to improve predictability such as seeking an optimal reliability cutoff point or excluding patients with papilledema either did not produce the intended practical outcome or further restricted its applicability [21]. Querfurth et al. introduced simultaneous color doppler monitoring of the central retinal and ophthalmic arterial flow velocities to improve predictability. They reported a more significant correlation (R = 0.95, p < 0.005) with combined parameter approach [23]. Modified photoplethysmography provides a quantitative ICP prediction along a continuous scale. Fourier domain decomposition of the retinal vascular pulse amplitude enables a selective inclusion of pulse wave harmonics hence, increasing the signal-to-noise ratio in the data set. In a recent study, we described a physiological model to estimate ICP where a plot of the induced intraocular pressure (x-axis) against the retinal venous pulse (y-axis) was measured using Modified Photoplethysmography. Predicted ICP was plotted from the x-axis intersection extrapolated from the peak retinal venous pulsation amplitude. A mean absolute error of 3.0 mmHg was achieved using this technique [29]. However, the XGB approach adds unique advantages: 1) The machine learning decision tree model addresses the heteroscedastic data structure directly and allows the model fit to be optimized as cases are added to the dataset thereby futureproofing the model's performance. 2) An improved predictive analysis by a reduction in the mean bias to 0.034±1.80 cm water (0.025±1.32 mmHg). 3) The ability to generate an independent prediction from the retinal arteries and veins and from all points from the retinal vessels in the image field enables visualization of the distribution of the predictions in the form of a peak density plot and to draw comparisons from the independent outputs. Therefore the XGB approach provides a method of internally validating the prediction from each case. However, papilledema may reduce the predictive accuracy when the renormalization of the ICP is followed by a lag in regression of the optic nerve changes [87], this may impact diagnostic accuracy of posttreatment serial Modified Photoplethysmography imaging tests.
Artificial intelligence classification methods that use clinical, electrophysiologic, or radiographic data to discriminate between normal and high ICP have been described [88][89][90]. Neural network classification models have yielded a total accuracy ranging between 70.2±4.5% to 92.05±2.25%. Golzan et al. estimated ICP using a neural network regression model derived from retinal venous pulse amplitude measured using the Dynamic Vessel Analyser, which is the only commercially available device used to quantitatively measure the retinal vascular pulse parameters. This device uses arbitrary units rather than frequency domain decomposition of the pulse wave. They reported a mean square error of ±1.27 mmHg in ICP prediction. However, they did not validate the model predictions on external cases, therefore, only internal validation results were reported. Moreover, only the venous pulse amplitude was used in the model prediction the arterial pulse characteristics was not taken into consideration in the model. In contrast to our method, the pulse amplitude was quantified using empirical units rather than a frequency domain decomposition strategy [91]. In our study, the XGB approach provided unique advantages over other machine learning options, hyperparameter tuning using Bayesian optimization (Table 1), mitigated model overfitting, this is the case where model output is not generalizable due to close fit to the sample data, a known limitation in decision tree machine learning algorithms. The wide normal ICP range (2-25 cm water) [39,92] and the predominance of non-linear dynamics in the interactions between ICP, IOP, and the retinal vascular pulse [28] may both contribute to data heteroscedasticity, which is the non-constant distribution of the error term of the predictors. This was addressed by the unique ability of a decision tree regression to generate a prediction without specifying an average structure for the model. The specific pruning algorithm is a significant factor in addressing the heteroscedasticity [93]. Moreover, the XGB analysis approach requires minimal data preprocessing and neither normalization nor scaling are necessary [94,95].
Intracranial pressure predictions derived from the venous model demonstrated better agreement with measured ICP compared to that provided by the arterial model. This may be due to the higher amplitude pulsation in the retinal veins compared to the arteries and the difference in venous pulsation amplitude between the ICP n and ICP h groups. Functional and structural differences between the arterial and venous systems, which in turn result in differences in wall tension and compliance may have played a role. Blood vessel walls consist mainly of water (70%), which is inelastic and incompressible, the remaining structure consists of a mesh of fibers with elastic properties. The fibrillary material consists of collagen, elastin, and smooth muscle cell layer in various percentages depending on the type and location of the vessel [96,97]. Heterogeneity of structural building blocks of the vessel wall results in a radiustension curve with non-linear characteristics. Whereas collagen fibers demonstrate a steep length-tension relationship and are closer to the linear relationship as expected from Hooke's law, elastin has a flat length-tension relationship [98], and that of smooth muscle is modified by the level of contraction [96]. Additionally, there are geometric factors that augment the venous response to changes in transmural pressure, at physiological ranges, a relatively small increase of transmural pressure from 0 to 10 mm Hg increases the venous volume by � 200%, which reflects a change in geometry from ellipsoidal to circular associated with an increased cross-sectional area. At supra-physiological venous transmural pressure (>40mmHg) there is a minimal increase in the relative volume as a result of the increase in compliance. In contrast, the arterial wall shows a smaller slope of the gradient between relative volume and transmural pressure accounting for the higher tolerance to transmural pressure [98]. Therefore, the arterial wall has curvilinear compliance over a range of transmural pressures, whereas venous walls have a semi-sigmoidal range of compliance. These structural and consequently functional differences may underlie the differences in vascular response to raised ICP and consequently may impact model predictive accuracy. Models derived from each vascular system misclassified one of the hold-out test cases with high ICP (case-2 in the arterial model and case-7 for the venous model). Additionally, both models predicted raised ICP for case 1 when the measured ICP was normal. The agreement between the predictions (�25 cm water) derived from the XGB models from both vascular systems in case-1 (Table 5) raises the possibility that the predicted ICP may represent the actual ICP. This scenario highlights the challenges of lumbar puncture manometry. Studies based on continuous monitoring report that ICP generally measures between 5 -10 mm Hg above atmospheric pressure [99]. However, this value may fluctuate within a range of 30 cm water over 12 hours [100]. Therefore, a time delay between the lumbar puncture and the ophthalmic assessment may have resulted in a spurious outcome in this case. Furthermore, a precise ICP estimate requires a consistent manometry technique where the needle entry point needs to be on the same level as the midline of the spine, which should also be at the same level as the patient's head. To complicate matters leakage around the needle can result in an erroneous estimate [6]. Further consideration needs to be given to the fact that ICP is estimated and measured from different levels. A historic assumption intracranially measured ICP (EVD-ICP) is equal to opening pressure measured at the lumbar spine (LP-ICP) [101,102]. In recent years, this matter has been a subject of debate where some studies have documented agreement between intracranial and lumbar measurements [103], correlated [104,105], and others have disputed both agreement and correlation [106,107]. Among these studies, Lenfeldt et al. experimentally raised the ICP using computerized lumbar infusion in subjects with communicating hydrocephalus. Simultaneously measured LP-ICP and Brain tissue ICP demonstrated a strong correlation (R 2 = 0.98) between the two ICP measurements throughout a pressure range of 0-600 mm water, the mean±standard deviation was -10±29 mm water. Therefore, LP-ICP correlates strongly with brain tissue-ICP in the absence of pathological obstruction to the cerebrospinal fluid flow [105].
Both feature importance (Figs 7 and 8) and SHAP summary plots (Figs 9 and 10) demonstrated that IOP i was the most significant parameter in both the arterial and venous models. Moreover, there was a positive correlation between IOP i and the impact on both XGB models. This indicates that ophthalmodynamometry is an important component of our image acquisition methodology. Moreover, the findings are consistent with the high correlation between retinal venous opening pressure and ICP described in ophthalmodynamometric studies in human [20-23] and animal studies [18,19,[108][109][110]. Changes in the venous system from raised ICP arises from elevated downstream cerebral venous pressures and external pressure on the subarachnoid segment of the central retinal vein, which increases vascular resistance [111]. The literature has been less clear on the changes in retinal arterial circulation in this context. Moreover, it is intriguing that arterial pulsation parameters provided information sufficient to generate a prediction in spite of the less clinically recognised changes in the retinal arterial system with raised ICP. Our study demonstrated an increase in the arterial wall pulsation amplitude in the ICP h group albeit at a narrower range than that of the venous system (Fig 4). In a recently published study, in addition to the retinal vascular pulse amplitude antipodal effect in response to abnormally elevated intracranial pressure, where a decrease in retinal venous pulse was accompanied by an increase in the arterial pulse proportionate to raised ICP [28]. We postulated that this effect may represent functional differences in the transmural pressure gradient between the retinal arteries and the retinal veins, where the former exhibits a higher range than the latter, therefore the retinal artery demonstrates pulsations at higher induced IOP i values, this combined with raised venous pressure and reduced venous compliance and possibly augmentation of the retinal arterial pressure wave in association with elevated ICP. Other compensatory mechanisms known to modify the cerebrovascular dynamics through auto-regulatory control have been described in the retinal arterial system as well [112][113][114][115]. However, the role of these mechanisms in our observations remain unresolved. We found that the SHAP values suggested that ocular laterality was a prominent factor in the predictive model, particularly that derived from the arterial system, it could be hypothesized that hemodynamic phenomena depend on physiologically inherent asymmetry in vascular dynamics play an important role in the optimization of cardiovascular functions [116,117]. Recent studies have demonstrated the association between ocular laterality on retinal vascular occlusions [118,119]. Although there are anatomic reasons for retinal arterial occlusions pertaining to embolic etiologies, the reason why asymmetry plays a role in the arterial rather than the venous system remains unclear.
A major limitation of machine learning algorithms is the black-box problem, given that the nature of machine learning is based on accuracy-driven performance metrics, these models will likely continue to become even more opaque in the future, especially machine learninggenerated ensembles of decision trees [120][121][122][123]. More recently, software packages like DALEX [124], breakdown [125] and XGBExplainer [126] have made some gains in terms of XGB model interpretability. A fundamental restriction inherent in all tree-based models is the inability to extrapolate target values beyond the range of the training data when making predictions, this is in contrast to linear models that can extrapolate predictions beyond the range of the training dataset [127]. Our study included a small number of participants, a larger dataset would likely offer improved generalizability. Image acquisition requires subject cooperation, which hinders the technique's applicability in patients with cognitive or significant neurological deficits, therefore future developments would include a compact imaging system.

Conclusion
The venous XGB model showed higher predictive accuracy compared to the arterial. Among the nine evaluated features IOP i was the most important model feature in both vascular systems. Although the venous median predicted ICP showed the highest agreement with measured ICP using lumbar puncture, an independent prediction derived from the retinal arteries can provide a system of internal validation for the result. Previous studies have not considered the retinal arterial system for this purpose. The prediction model can potentially be integrated into a neurological clinical decision algorithm to evaluate the indication for lumbar puncture.